'''
Created on 2010-11-4

@author: yun_hua

Reference: 1. Numerical Heat Transfer, 2rd, Wenquan Tao, Xi'an Jiaotong University
           2. Elements of Chemical Engineering, 4th, H.Scott Fogler
           
           Using Finite Volume Method(FVM), and the first-order upwind(FUW) method.
           
           In the mass balance equations, the component A's variation is represented by rhoA / rhoAll.
           So the general diffusion coefficient is diffusivity * rhoAll (here we assume rhoAll is constant).
           
           About Simpson integration:
           Suppose xk = a + k*h, k = 0,1,...,2M, [a, b] is divided into 2M spaces while h = (b - a)/(2*M)
           S(f, h) = h / 3 * (f(a) + f(b)) + 2 * h / 3 * Sum(f(x2k))k=1 to M -1 + 4 * h / 3 * Sum(f(x2k-1))k=1 to M
'''

import math
import numpy as N
import phyconst

X_INDEX = 0
Y_INDEX = 1
Z_INDEX = 2

class TubularReactor(object):
    """
    The tubular reactor model simulation using FVM method.
    """
    def __init__(self, grid, method, param):
        """
        Constructor method.
        grid: the grid for calculation.
        method: the method used to express the convection discretion: Upwind, QUICK 
        param: the parameter factory
        """
        if grid is None or grid.size == 0:
            raise ValueError, "The grid is empty or null."
        
        # Create the solution array for the temperature and the mass density percentage
        self.grid = grid      
        self.radius = self.grid[-1, -1, 1]
        self.len = self.grid[-1, -1, 0]
        self.tarray = N.zeros((self.grid.shape[0], self.grid.shape[1])) # The temperature array for the calculation
        self.marray = N.zeros((self.grid.shape[0], self.grid.shape[1])) # The (m)ass percentage array for the calculation
        
        # Convert the material properties to an array
        prop = param.get_componentprop()
        self.prop = N.zeros((len(prop.keys()), len(prop.values()[0])))
        self.prop[0, :] = prop['A-PropyleneOxide']
        self.prop[1, :] = prop['B-Water']
        self.prop[2, :] = prop['M-Methanol']
        self.prop[3, :] = prop['C-PropyleneGlycol']   
        """ FVM parameters
             big deltax
            |         |
         W--|----P----|--E
            |         |
         |deltax |
        1. For the uniform grid, bigdeltax == deltax
        2. For the cylinderical coordinates, deltay = delta_r * radius
        """ 
        self.deltax = abs(self.grid[0, 1, 0] - self.grid[0, 2, 0])
        self.bigdeltax = self.deltax
        self.deltay = abs(self.grid[1, 0, 1] - self.grid[2, 0, 1])
        self.bigdeltay = self.deltay
        self.deltaCp = (self.prop[3, param.cpIndex] * self.prop[3, param.stioIndex] - \
                    self.prop[1, param.cpIndex] * self.prop[1, param.stioIndex] - \
                    self.prop[0, param.cpIndex] * self.prop[0, param.stioIndex]) / self.prop[0, param.stioIndex]
        self.deltaHrx = (self.prop[3, param.hxIndex] * self.prop[3, param.stioIndex] - \
                    self.prop[1, param.hxIndex] * self.prop[1, param.stioIndex] - \
                    self.prop[0, param.hxIndex] * self.prop[0, param.stioIndex]) / self.prop[0, param.stioIndex]                     
        # Other initial values
        self.actenergy = param.get_actenergy()
        self.t0 = param.get_t0()
        self.t0base = param.get_t0base()
        self.t0wall = param.get_t0wall()
        self.k0 = param.get_k0()
        self.ke = param.get_ke()
        self.uht = param.get_uht()
        self.de = param.get_de()
        
        self.f0 = sum(self.prop[:, param.volumeRateIndex])
        self.ux = self.f0 / (phyconst.pi * self.radius ** 2)
        self.cpm = sum(self.prop[:-1, param.moleRateIndex] * self.prop[:-1, param.cpIndex]) / self.f0
        self.rho = sum(self.prop[:-1, param.moleRateIndex] * self.prop[:-1, param.moleWeightIndex]) / self.f0
        self.mA0 = self.prop[0, param.moleRateIndex] * self.prop[0, param.moleWeightIndex] / \
                         sum(self.prop[:-1, param.moleRateIndex] * self.prop[:-1, param.moleWeightIndex]) # percentage: component A mass / the total mass 
        
        # Set the initial condition including the boundary.
        self.init()
        
    def getcoffarray(self, lp, mp):
        """
        Get the coefficient array for the FVM calculation.
        lp: The l-direction position
        mp: The m-direction position
        """
        return N.zeros(1)
    
    def init(self):
        """
        Set the initial conditions including boundaries, temperature and the component mass percentage.
        """
        self.tarray += self.t0
        self.marray += self.mA0
        self.wallboundypos = self.grid.shape[0] - 2
        self.axisboundypos = 1
        self.inletboundxpos = 1
        self.outletboundxpos = self.grid.shape[1] - 2
    
    def getboundary(self):
        pass
    
    def getmarray(self):
        """
        Get mass array.
        """
        return self.marray.copy()
    
    def setmarray(self, marray):
        """
        Set mass array.
        """
        self.marray = marray
    
    def gettarray(self):
        """
        Get the temperature array.
        """
        return self.tarray.copy()
    
    def settarray(self, tarray):
        """
        Set the temperature array.
        """
        self.tarray = tarray
        
    def updatetemp(self, xpos, ypos, t):
        """
        Update the temperature in the given position
        """
        self.tarray[ypos, xpos] = t
    
    def updatemass(self, xpos, ypos, m):
        """
        Update the mass in the given position
        """
        self.marray[ypos, xpos] = m
    
    def getgrid(self):
        """
        Get the node grid.
        """
        return self.grid.copy()
    
    def getaE(self, xpos, ypos, isheatflux):
        """
        Get the FVM method's coefficient for East point, using the first-order upwind scheme(FUS).
        xpos: the x-axis position
        ypos: the y-axis position
        isheatflux: true if it is calculating for heat flux, otherwise the mass flux
        return: the East point factor
        """ 
        if isheatflux:
            diff = self.ke
        else:
            diff = self.de * self.rho
        if not xpos == self.outletboundxpos:
            return (diff / self.deltax) * self.bigdeltay * self.grid[ypos, xpos, Y_INDEX] 
        else:
            return 0.
    
    def getaW(self, xpos, ypos, isheatflux):
        """
        Get the FVM method's coefficient for West point, using the first-order upwind scheme(FUS).
        xpos: the x-axis position
        ypos: the y-axis position
        isheatflux: true if it is calculating for heat flux, otherwise the mass flux
        return: the West point factor
        """
        if isheatflux:
            diff = self.ke
            fluxrho = self.cpm * self.ux
        else:
            diff = self.de * self.rho
            fluxrho = self.rho * self.ux
        return (diff / self.deltax + fluxrho) * self.bigdeltay * self.grid[ypos, xpos, Y_INDEX]
        
    def getaN(self, xpos, ypos, isheatflux):
        """
        Get the FVM method's coefficient for North point, using the first-order upwind scheme(FUS).
        xpos: the x-axis position
        ypos: the y-axis position
        isheatflux: true if it is calculating for heat flux, otherwise the mass flux
        return: the North point factor
        """
        if isheatflux:
            diff = self.ke
            fluxrho = 0. # The Y(radius)-direction flux is omitted.
        else:
            diff = self.de * self.rho
            fluxrho = 0. # The Y(radius)-direction flux is omitted.
        if not ypos == self.wallboundypos:
            return (diff / self.deltay + fluxrho) * self.bigdeltax * self.grid[ypos, xpos, Y_INDEX]
        else:
            return 0.
    
    def getaS(self, xpos, ypos, isheatflux):
        """
        Get the FVM method's coefficient for South point, using the first-order upwind scheme(FUS).
        xpos: the x-axis position
        ypos: the y-axis position
        isheatflux: true if it is calculating for heat flux, otherwise the mass flux
        return: the South point factor
        """
        if isheatflux:
            diff = self.ke
            fluxrho = 0. # The Y(radius)-direction flux is omitted.
        else:
            diff = self.de * self.rho
            fluxrho = 0. 
        if not ypos == self.axisboundypos:  
            return (diff / self.deltay + fluxrho) * self.bigdeltax * self.grid[ypos, xpos, Y_INDEX]
        else:
            return 0.
    
    def getaP(self, xpos, ypos, aE, aW, aN, aS, isheatflux):
        """
        Get the FVM method's coefficient for the Center point, using the first-order upwind scheme(FUS).
        xpos: the x-axis position
        ypos: the y-axis position
        valuearray: the mass or temperature value array
        boundarytype: type 0, not boundary; type 1, Dirichlet condition; type 2, Neumann condition.
        boundaryflux: the boundary mass or heat flux density.
        boundaryvalue: the boundary value.
        isheatflux: true if it is calculating for heat flux, otherwise the mass flux
        """
        vunit = 0.
        sp = 0.
        if ypos == self.wallboundypos and isheatflux:
            return aE + aW + aN + aS + \
                (self.bigdeltax * self.grid[ypos, xpos, Y_INDEX] / (1. / self.uht + self.deltay / self.ke)) \
                - sp * vunit
        else:
            return aE + aW + aN + aS - sp * vunit
    
    def getsource(self, xpos, ypos, isheatflux):
        """
        Get the FVM method's coefficient for the source, using the first-order upwind scheme(FUS).
        xpos: the x-axis position
        ypos: the y-axis position
        isheatflux: true if it is calculating for heat flux, otherwise the mass flux
        """
        vunit = self.bigdeltax * self.bigdeltay * self.grid[ypos, xpos, Y_INDEX]
        sc = -self.k0 * \
                math.exp(self.actenergy / phyconst.Rgas  * (1. / self.t0base -  1. / self.tarray[ypos, xpos])) * \
                self.marray[ypos, xpos] * self.rho # rate law
        if isheatflux:
            sc /= self.prop[0, 0] # Convert the density to the concentration: kg/m**3 -> mol/m**3 
            sc *= self.deltaHrx + self.deltaCp * (self.tarray[ypos, xpos] - self.t0base)
        # Tubular's wall boundary condition for the heat transfer
        if isheatflux and ypos == self.wallboundypos:
            sc += self.bigdeltax * self.grid[ypos, xpos, Y_INDEX] * self.t0wall / vunit / \
                (1. / self.uht + self.deltay / self.ke)
        # The wall boundary for the mass, the uht = 0
        
        return sc * vunit
        
    
    def get_deltax(self):
        """
        Get small delta x.
        """
        return self.deltax
    
    def get_deltay(self):
        """
        Get small delta y.
        """
        return self.deltay
    
    def fillboundary(self):
        """
        Fill the boundary values according to the boundary conditions.
        """
        self.marray[1:-1, -1] = self.marray[1:-1, -2] # The outlet boundary
        self.marray[0, 1:] = self.marray[1, 1:] # The axis boundary
        self.marray[-1, 1:] = self.marray[-2, 1:] # The wall boundary
        
        self.tarray[1:-1, -1] = self.tarray[1:-1, -2] # The outlet boundary
        self.tarray[0, 1:] = self.tarray[1, 1:] # The axis boundary
        self.tarray[-1, 1:] = self.tarray[-2, 1:] # The wall boundary
        
        q = (self.t0wall - self.tarray[-2, 1:]) / (self.ke / self.uht + 0.5 * self.deltay / self.ke)
        self.tarray[-1, 1:] = self.t0wall - q * (self.ke / self.uht)
        
    def getconvprofile(self):
        """
        Get the mean conversion profile along the the x-axis, using the Simpson integration method.
        """
        h = self.deltay
        conv = []
        radiusgrid = self.grid.shape[0] - 2
        rarray = (N.arange(0., radiusgrid + 1) / radiusgrid) * self.radius
        for nx in range(self.grid.shape[1]):
            marray = (self.marray[1:-2, nx] + self.marray[2:-1, nx]) / 2.
            integ = (self.rho / self.prop[0, 0]) # Component A's density to mole concentration
            integ *= self.ux
            integ *= h / 3. * \
                    (self.marray[0, nx] * 0. + self.marray[-1, nx] * self.radius) + \
                    2. * h / 3. * (marray[1::2] * rarray[2:-1:2]).sum() + 4. * h / 3. * (marray[0::2] * rarray[1:-1:2]).sum()
            conv.append(1. - 2. * phyconst.pi * integ / self.prop[0, 6])
        return conv
    
    def getconcentrationprofile(self):
        """
        Get the mean concentration profile along the x-axis, using the Simpson integration method.
        """
        h = self.deltay
        con = []
        radiusgrid = self.grid.shape[0] - 2
        rarray = (N.arange(0., radiusgrid + 1) / radiusgrid) * self.radius
        for nx in range(self.grid.shape[1]):
            marray = (self.marray[1:-2, nx] + self.marray[2:-1, nx]) / 2.
            integ = (self.rho / self.prop[0, 0]) # Component A's density to mole concentration
            integ *= self.ux
            integ *= h / 3. * \
                    (self.marray[0, nx] * 0. + self.marray[-1, nx] * self.radius) + \
                    2. * h / 3. * (marray[1::2] * rarray[2:-1:2]).sum() + 4. * h / 3. * (marray[0::2] * rarray[1:-1:2]).sum()
            con.append(2. * phyconst.pi * integ / math.pi / self.radius**2 / self.ux)
        return con
        